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Abstract 

Background: The concept of disaster surge has arisen in recent years to describe 
the phenomenon of severely increased demands on healthcare systems resulting 
from catastrophic mass casualty events (MCEs) such as natural disasters and terrorist 
attacks. The major challenge in dealing with a disaster surge is the efficient triage 
and utilization of the healthcare resources appropriate to the magnitude and 
character of the affected population in terms of its demographics and the types of 
injuries that have been sustained. 

Results: In this paper a deterministic population kinetics model is used to predict 
the effect of the availability of a pediatric trauma center (PTC) upon the response to 
an arbitrary disaster surge as a function of the rates of pediatric patients' admission 
to adult and pediatric centers and the corresponding discharge rates of these 
centers. We find that adding a hypothetical pediatric trauma center to the response 
documented in an historical example (the Israeli Defense Forces field hospital that 
responded to the Haiti earthquake of 2010) would have allowed for a significant 
increase in the overall rate of admission of the pediatric surge cohort. This would 
have reduced the time to treatment in this example by approximately half. The time 
needed to completely treat all children affected by the disaster would have 
decreased by slightly more than a third, with the caveat that the PTC would have to 
have been approximately as fast as the adult center in discharging its patients. Lastly, 
if disaster death rates from other events reported in the literature are included in the 
model, availability of a PTC would result in a relative mortality risk reduction of 37%. 

Conclusions: Our model provides a mathematical justification for aggressive 
inclusion of PTCs in planning for disasters by public health agencies. 



Background 

In the modern era, humanity has spread across and settled all habitable areas of the 
globe, thereby greatly increasing potential exposures to catastrophic events, whether 
natural or manmade, as demonstrated most recently by the 2010 Haiti earthquake [1] 
as well as the tragic earthquake, tsunami and nuclear disaster that devastated Japan in 
March, 2011 [2]. It is imperative that planning be undertaken to deal effectively with 
the vast number of injured survivors. These conditions can be described as a disaster 
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surge, which can be thought of as an unusually high fluctuation over and above the 
normal background rate of patient utilization of medical services [3-12]. Multiple 
strategies have been proposed to maximize patient throughput and efficiency of 
resource utilization under surge conditions, and the overall consensus is that detailed 
planning for various disaster contingencies is the key to this process. 

Because of the random, stochastic nature of disaster events, this planning can be 
greatly aided by simulation. A considerable amount of work has been done in modeling 
disaster surges and the response of health systems to them [13]. More generally, a 
patient population having to wait for medical triage and treatment can be thought of as 
a problem in queueing theory [14-17]. This field grew out of A. K. Erlang's pioneering 
approach to modeling demand for telephone service in the early 20 th century [18,19], 
and has been applied to a diverse range of problems including not only telecommunica- 
tions, but airport and automobile traffic patterns, other service industries, and hospital 
and factory design [20-22]. If the length of the queue is long, then its behavior can often 
be approximated to that of a continuous variable, thereby simplifying the mathematics 
greatly. This approach results in what are referred to in the queueing theory literature as 
fluid models [23-25], and can be used for predicting the behavior of, for example, queues 
for service from a call-in center [26] . It has also been shown that if a system satisfies the 
Markov property, that is, if its future behavior depends only on its current state, then its 
behavior can be approximated deterministically by simple ordinary differential equations 
(ODE's) [27,28]. While more complicated stochastic methodologies such as Monte 
Carlo simulation have been successfully used in modeling the response to a patient 
surge [29,30], the simplicity of the ODE approach has motivated the use of kinetic or 
compartmental models for such problems [31]. In this method, the population evolves 
from an initial state to a number of subsequent states with each state change having a 
rate constant. This approach has also long been used in physics and chemistry to model 
reactions and series of reactions, as well as in population biology [32-34]. Here, we make 
use of this mathematically elementary and well-established approach to predict the 
behavior of pediatric and adult populations after a mass casualty event, with and without 
the availability of a facility specifically designed to treat children. 

A significant proportion of disaster victims are children, who have unique physiology, 
patterns of injury, and psychosocial needs in such settings [35]. Studies have shown that 
the availability of a pediatric trauma center (PTC) would probably improve the overall 
response to a mass casualty incident, but the available data are sparse [36]. In the absence 
of more extensive data, in this paper we use a population kinetics approach to estimate 
the effect of the availability of a pediatric trauma center upon the rates of admission and 
discharge of a disaster surge population by extrapolating from historical data. We find 
that the initial rate of discharging patients from the PTC early in the surge is the dominant 
influence on the time needed to fill the hospitals maximum bed capacity as well as on the 
time needed to definitively treat and discharge all patients in the surge. On the other 
hand, the PTC admission rate and the rate of discharging patients once the PTC is full are 
the most important factors in determining the time needed to admit the entire surge. We 
then add historical mortality rates to our model and calculate the reduction in deaths that 
would be conferred by a PTC. We conclude that within the limits of our model, the avail- 
ability of a PTC would greatly enhance the response to a disaster as measured by the total 
time needed to appropriately triage and treat the surge population. 
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Methods 

I. Approach 

Before describing the details of our model, we shall first solve a simpler problem that 
will provide its mathematical underpinnings. We begin by assuming that an unspeci- 
fied disaster instantaneously produces an initial surge population. This scenario is a 
good approximation for a subset of mass casualty events (MCEs) that occur suddenly 
without appreciable buildup or exposure time, such as bombings, earthquakes, or air- 
plane crashes. (The more general case, where there is a delay between the inciting 
event and the onset of the surge, is mathematically more complicated, requires more 
unknown parameters than the current scenario, and is developed for completeness in 
Appendix A.) This population, which we shall denote by N s (t), is defined at time zero 
to be N s (t = 0) = N 0 , and changes as it is admitted to a trauma center into a popula- 
tion N a (t) of admitted patients with rate k a , which in turn can become a population of 
N d (t) discharged patients with rate k^: 

N s hNA N d (1) 

Appropriate estimates for /c^and /c^will be discussed later when we apply our model 
to real-world historical data. We note that "discharge" would include mortality in this 
scheme, as no explicit provision is made for categories of discharge (discharged to 
home, discharged to a long term care facility, deceased, etc.). Equation 1 governs the 
behavior of the surge population as patients transition to being admitted and treated, 
and ultimately discharged; this behavior is described mathematically by a set of three 
coupled first-order differential equations: 

-r- = ~k a N s (2) 
at 

^ = k a N s -k d N a (3) 
at 

IT = kdNa (4) 

To solve Eqs. 2-4, we require the boundary conditions: 

N s (t = 0) = N 0 (5) 

N s (t -> oo) = 0 (6) 

N a (t = 0) = N a {t -> oo) = 0 (7) 

N d {t = 0) = 0 (8) 

N d (t -> oo) = No (9) 

Eqs. 5 and 6 state that the number of surge patients begins at N 0 , and decays to zero 
at long times since all patients are admitted and discharged. Eq. 7 reflects the fact that 
there are no patients admitted at time zero, and at long times all admitted patients 
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have been discharged. Eqs. 8 and 9 therefore state that there are no discharged patients 
at time zero, while at long times the entire population has been discharged. 

We can now solve the system of equations 2-4. Equation 2 can be solved by direct 
integration, and applying the boundary conditions 5 and 6 gives: 

N s (t) =N 0 e- kat (10) 

Eq. 10 can be substituted into equation 3, yielding with some rearrangement: 

< ^L+k d N a =k a N 0 e- k ° t (11) 
at 

Multiplying Eq. 11 by the integrating factor exp(/c^), integration and application of 
boundary condition (7) gives: 

^a[t) = N 0 -\- (e- k ^ - e~ k A (12) 

This can be substituted into Eq. 4, which after direct integration and application of 
boundary conditions (8) and (9) gives 

N ^ = N ° ( v\ e ~ kdt ~ v\ e ~ kat + 0 < 13 > 

\kd-ka k d -k a J 

Figure 1 shows a schematic of the behavior of the populations N s , N ai and A/^as 
described by Equations 10, 12 and 13. No units are shown here for the sake of concep- 
tual clarity; quantitative results are shown in the Results section. The surge population 
decays with typical single exponential behavior; the admitted population rises to a 
maximum and decays, and the discharged population exhibits an exponential rise. 

II. Maximum Capacity Model 

At this point, we note that the model as currently formulated has a limitation in that 
no provision is made for the maximum capacity of the trauma center. In other words, 
the maximum value of N a (t) predicted by Eq. 12 is a function only of N 0 , /c^and k& 
with no dependence on the number of available beds in the center. To see this, N a (t) 
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can be maximized by setting its derivative equal to zero and solving this expression for 
U which gives 

t Nmax = 1 In ( ^ ] (14) 



h-h \k a 

This value for t is inserted back into Eq. 12, giving 

-K -kd 



N " ax = N o i -^- 
kd - K 



h)kd-k a _(*d\k d -k a 

K) \k a ) 



(15) 



which is a function only of N 0t /c^and kj, and has no relation to any real-world hospi- 
tal bed capacity. 

This limitation can be overcome by modifying the model with some intuitive assump- 
tions. First, we shall identify our trauma center's intrinsic maximum capacity as A/^ max , 
and assume that the surge population behaves just as we have described above until 
A^reaches A/^ max . This maximum census is not equal to the total number of beds in the 
trauma center, but rather its surge capacity over and above normal operations, or equiva- 
lently the fraction of its beds allotted in the center's planning for an MCE [37]. After 
N a max is reached, we assume that the center will remain at maximum capacity until the 
surge is exhausted. That implies that the admission and discharge rates are equal during 
this period. Next, we assume that the admission rate will be somewhat lower after the 
trauma center is full compared with early times, as during this period many of its surge 
beds will be occupied with critically injured patients. Finally, we assume that once 100% of 
the surge has been admitted, the trauma center's discharge rate will return to that prior to 
maximum patient load. We will call this modified model the "maximum capacity model." 
To formulate this modification of the model mathematically, it is helpful to define two 
times ti and t 2 as illustrated in Figure 2. At t lf the trauma center has reached its maximum 
capacity and can only admit a patient if another is discharged, i.e., t\ = t^ a max • This situa- 
tion persists until t 2) at which time the surge population has declined to zero and the 
trauma center can again discharge patients at the pre-MCE rate. In the language of 







m 





























*1 h t 

Figure 2 Behavior of surge (blue), admitted (red), and discharged (green) populations with time in 
the maximum capacity model described by Equations 16-18. Curves are not normalized or scaled. 
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queueing theory, t 2 is the time at which the queue has vanished. The resulting population 
behavior is shown in Figure 2, where again, no units are shown for conceptual clarity; 
quantitative data are shown below in Results. 

We are now ready to solve the necessary differential equations for the maximum 
capacity model. Overall, we have three separate regions in time with different behavior, 
defined by 

-k a N s 

k a N s - k d N a (16) 

k d N a 



0 < t < h 



dNs_ 
dt 
dN a 

dt 
dN d 

dt 



h < t < t 2 



dN s 

dt 
dN a 

dt 
dN d 

dt 



= -k f 



(17) 



dN s 

t>t 2 : - = 0 

^ = -hN a (18) 
at 

dN d U 

—j— = hN a 
dt 

All variables and parameters have the same meanings as previously defined, except 
for a single new parameter k' that describes the admission and discharge rates during 
the period of time between t x and t 2 when the trauma center is operating at maximum 
capacity. We again note that k' will likely be less than either k a ov kj, as both admis- 
sions and discharges will be slower once the trauma center is filled with critically 
injured surge patients. 

The solution of Eq. 16 is identical to Equations 10, 12 and 13. However, at t x the 
system's behavior changes to conform to Equation 17, giving 

h < t < t 2 : N s = N sA - k'(t - h) 

N a = N aA (19) 

N d = N ditx +k'{t-t l ) 

where N a , tl = A/^ max . At t 2 , the entire surge population has been exhausted, and the 
system's behavior changes to that entailed by Eq. 18, the solution of which is 

t > t 2 : N s = 0 

Na=N a>t2 e- k ^ (20) 
N d =N dit2 +Na lt2 (l-e- k ^- t ^ 
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We note that t x has already been determined to be equal to tMa max as defined in Eq. 
14 and again N a ,t 2 = N a ,ti = N a m ^. We are now also in a position to determine t 2 , 
which is the time when N a = 0. Eq. 19 then gives 

t 2 = h + ^ (21) 



III. Maximum Capacity with Pediatric Trauma Center Model 

Armed with Equations 16-18 we can now include the effect of an available pediatric 
trauma center in the maximum load model. We shall call what follows the "maximum 
capacity with pediatric trauma center model" We now assert that the initial N 0 disas- 
ter victims are composed of A 0 adults and P 0 pediatric patients, viz: 

N 0 = A 0 + P 0 (22) 

We also note that the total number of surge patients as a function of time is equal to 
the sum of the adult and pediatric subpopulations: 

N s (t) = P s (t) + A s (t) (23) 

where A 5 and P 5 now indicate the adult and pediatric cohorts of the surge, respectively. 
We then assume that adult patients are only admitted to adult trauma centers, while 
pediatric patients may be triaged and admitted to either adult or pediatric trauma cen- 
ters (PTCs); this assumption is similar to the approach taken by Perry and Whit in 
modeling call center capacity overloads [26], except that our case is asymmetric: adults 
are never triaged to PTCs in this model. These assumptions result in the following 
kinetic scheme: 

As hA a H A d (24) 

Ps ^ Paa ^ Pa ( 25 ) 

Ps ^ P ap H P d (26) 

where /c^and /^represent the rates of adult admission to and discharge from an 
adult center, k paa s.nd k pc i a the rates of pediatric admission to and discharge from the 
adult center, and /c^^and /c^the rates of pediatric admission to and discharge from 
the PTC. Similarly, A a (t) and A d (t) are the populations of admitted and discharged 
adults, while P aa (t) an d P a p{t) are the pediatric populations admitted to adult and 
pediatric centers, respectively, and P d {t) represents the discharged pediatric population, 
irrespective of the center at which they were treated. 

The differential equations entailed by Equation 24, boundary conditions, and their 
solution are identical to Equations 1-4 except for subscripts: 

A s (t) =A 0 e~ kaat (27) 
Aa[t) = A 0 - K \ (e-^ - e~ k A (28) 
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Kaa &ad ^aa / 

On the other hand, the coupled differential equations resulting from Equations 25-26 
are slightly different: 



dP s 
dt 



= -{kpaa + k pap )P s (30) 



dPaa 

^ = kp aa P s — kpd a P aa (31) 

dP 

-^=k pap P s -k pdp P ap (32) 

- kpdaPaa + kpdpPap (33) 

The differences arise from the fact that there are potentially different rates of admis- 
sion of pediatric patients to, and discharge of these patients from, the adult and pedia- 
tric trauma centers in the model. If the admission and discharge rates are equal, 
Equations 31-32 collapse into a single equation that is analogous to (1) and (24). The 
boundary conditions on (30) and (33) are the same as (5-6) and (8-9); those for (31) 
and (32) are identical to (7). At this point it is helpful to define: 

k = k paa + kp ap (34) 

Eq. 26 essentially defines an effective or total admission rate constant for pediatric 
patients in the model. The solution to (30-33), though slightly more complicated, is 
obtained via the same algorithm that led to (10-13) and is as follows: 

P 5 (t) = P oe - kt (35) 
Paa{t)=P 0J ^^(e- kt -e-^A (36) 

kpda - k\ / 

Pa P {t) = P 0 (e~ kt - e-W) (37) 



kpdp - k 

Pd(t) = P 0 



k P aa c -hdat + e -kp dp t 
kpda ~ k kp dp - k 

kpdakpaa _ kpdpkpap \ + ^ 



k{k - k pda ) k(k - k pdp ) 

Equations 35-38 describe the behavior of the pediatric cohort of the surge prior to 
the maximum load times for the adult and pediatric trauma centers. The behavior will 
change to one similar to Eq. 19 after these maxima are reached. However, there is no 
longer a single time for the maximum load, but rather separate ones for the adult and 
pediatric centers. Moreover, the behavior between the maximum load for the faster 
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facility and that of the slower one requires a separate set of differential equations to be 
solved. To show this, Figure 3 displays the qualitative behavior of the model we are 
about to derive, with the necessary boundary conditions and regions where each set of 
differential equations holds (I, II, III and IV) noted on the figure. The time for maxi- 
mum load for the center that admits and discharges patients more rapidly is analogous 
to Eq. 14. For the purposes of developing the model, we shall assume that the adult 
center is faster, but the derivation proceeds identically if the opposite assumption is 
made, except for the subscripts on the parameters; this issue is discussed further in 
additional files 1 and 2. We shall also omit the constant prefactor P 0 from all the equa- 
tions that follow, since it can be added back in after the derivation is complete with no 
loss of generality. Therefore: 
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The equations governing the system's behavior in region II of Figure 3 are 

dP s , 

h,a < t < t\ t p ! —j^- = —k a — kp a pP s 

dP 

aa = Q 

dP ap m 
^ = kpapPs ~ kpdpPap 

dP s , 

~j~ - k a + KpdpP a p 

The solution to these equations is 

h,a < t < t\ t p : 

kpap 
Paa = C 
ka 



1 kpdp V ) ( 41 ) 



P IW = ^-( g-**(t-tu) - l) 

P a =(h- Y + y~) (l - e-^- 1 ^) 

+y b± Ci_ e -W*-*uA +B 
kpap V / 

where the constants are given by 



a = A + 

Akp a p + k a 



kpap 



kpdp ~ kpap 



A = e~ ktl * 

r - kpaa ( A p-kpdatlA 

C ~l^k\ A - e ) (42) 

kpda ~ k kpdp ~ k 

kpdafcpaa kpdpkpap . ^ 



fe(fe - fepdfl) fe(fe - kpdp) 

R= J?P^( A _ e -k pdpKa \ 
kpdp ~k\ J 

The time for the maximum load D on the slower center (the pediatric trauma center 
in this derivation) can be found by maximizing the expression for P ap m (41). Thus the 
expressions for t lfP and D are: 

1 7 / Otkp a p y 

= I TT~ ln 



kpap — kpdp \fikpdp fi / 

D = P ap (K P ) = (e-^~^ - l) (43) 

+ Y ^e~kpap(h,p-ti,a) _ e -kpdp{h,p-h,a)\ 

+]-[e~kpdp {h.p-Ka) 
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with p given by: 



P = H — y + — — (44) 

Kpdp 

B and F are also obtained by substituting t ltP into the appropriate expressions in (41): 

k ' 

B = — + ae-kpapfap-Ka) 

kpap 

H-y + y-^j (l - e-Kpdpih.p-h.a)^ ( 45 ) 

+ yb±( 1 - e -hap{h,p-t ll a)\ +E 

kpap V ' 

After t h p , the slower pediatric center is also at its maximum capacity, and both cen- 
ters can admit a patient only if another is discharged. The kinetics then becomes zer- 
oth order similarly to Eqs. 17 and 19, so in region III of Figure 3, the governing 
equations are 

dP s 

h lP <t<t 2 : — = —k 
at 

dPaa = dPa L=o (46) 

dt dt 
dP d 
dt 

where n = k a ' + k p \ the sum of the discharge rates of the adult and pediatric centers 
during the period of time while they are at maximum capacity (region III). The solu- 
tion to (46) is 

h,p<t<t 2 : P s =B-K(t-t hp ) 

P aa — C 



Pap = D 



(47) 

Pd = F + K{t - t llP ) 

The value for t 2 , the time when the surge is exhausted, is just the ^-intercept of the 
line describing P s in region III. This can also be used to find G: 

t2 = t ^ + ? (48) 
G = F + K{t 2 - t hp ) = F + B 

Because the surge cohort has vanished at t 2 , there is no more external load upon 
either trauma center, so we again make the assumption that each center can resume 
discharging patients at the pre-MCE rate as in Eq. 18. This is only a simplifying 
assumption, as arbitrary rates could be assumed with no effect on the derivation except 
for subscripts. The differential equations for region IV are then 
dP s 

t>t 2 : — = P s = 0 
dt 

— — - -h , P 

dp a (49) 

~ir = ~ kpdpPap 

- kpdaPaa + kpdpPap 
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All boundary conditions are known and illustrated in Figure 3, including that at long 
times P^must be unity. Therefore the solution to (49) is: 

t > t 2 : P s = 0 

P aa = Ce- k P da{t - t2] 

i r , (50) 

P ap = De-kpdpit-ti) 

p d = l _ Q e -kpda(t-t 2 ) _ £) e -kpdp{t-t 2 ) 



Results 

I. Application of the maximum capacity model to an historical example 

Equations 40-50 now allow us to examine the behavior of the pediatric surge popula- 
tion P 0 under a variety of conditions. We begin by identifying the appropriate para- 
meters in the simpler maximum capacity model that define real-world timescales. We 
then proceed to work through an example of applying the model by considering litera- 
ture admission and discharge data from an historical disaster surge. We fit the equa- 
tions to these data, and then include the full maximum capacity with pediatric trauma 
center model to extrapolate the effect a pediatric trauma center would have had on the 
time necessary to treat the patients. 

There are several potentially observable parameters in the models presented here. 
The rates of admission and discharge in the initial and maximum capacity regimes are 
certainly observable in principle, but they are rarely reported as such. Also, the maxi- 
mum surge capacities NJ™**, C and D are available to disaster planners, but not usually 
reported directly. Rather, what is often available are the times of maximum load (t ± in 
the maximum capacity model, t ha and t ltP in the maximum capacity model with pedia- 
tric trauma center available) and the time at which the surge population has been 
completely dispositioned. The latter time does not correspond to t 2 , since the trauma 
centers are still full to capacity at this point. Rather, this is the time at which, in region 
IV of Figure 3, the discharged population has increased to very nearly unity. We note 
that it cannot be defined as the time that exactly 100% of the surge has been dis- 
charged, since the exponentials governing the behavior of the populations do not reach 
this value until infinity. Rather, we can define a time at which some specified fraction 
of discharges has been reached: we shall choose 99% and call this time t 99 . From Equa- 
tion 20, it follows that in the maximum capacity model t 99 is given by 

1 , / 0.99 -N dh \ , x 

t99 = t2 _ i n i d JL (51) 

kd V N a#t2 / 

However, in the maximum capacity with PTC model, Equation 50 is transcendental 
so t 99 cannot be solved for in closed form, but it can be found numerically. We note 
that our choice of the parameters t x and t 99 was motivated in large part by the avail- 
ability of such data in the literature, but also by the importance of t\ as a defining 
timescale of the behavior of populations in the model. On the other hand, we include 
t 2 primarily as a natural timescale of the model itself (where the surge or queue length 
vanishes and the system's deterministic behavior changes again) rather than as a 
descriptor of available historical data, and we examine the effect of varying it in the 
sensitivity analysis. Finally, the effect of including the explicit contribution of death 
rates for each population is derived in Appendix B. 
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The historical example we shall use in demonstrating the implementation of the model 
is that of an Israeli Defense Forces mobile field hospital that responded to the 2010 Haiti 
earthquake. In that case, 1111 patients were treated over a course of 10 days, and the hos- 
pital's maximum capacity of 60 to 72 beds was reached prior to 2 days of operation [1]. 
We chose this example because the disaster itself was sudden as required by our assump- 
tions, and because of the quality of the data available with which to fit our model in com- 
parison to other historical MCEs. We will begin by fitting the maximum capacity model 
(without a pediatric trauma center available) to these data. This amounts to solving a sys- 
tem of equations consisting of (13), (14) and (21) for k a , and k' with the historical data 
of ti = 2 days, t 99 = 10 days. Since we have three equations with two unknowns, the 
system is not uniquely determined and actually has two solutions, one for the case of k a >k- 
^and another for /c a </c d . To overcome this we must impose a constraint for t 2 : for this case 
we shall arbitrarily assume that it took approximately the same time to discharge all 
admitted patients once the surge was exhausted as it did for the hospital to reach maxi- 
mum capacity, that is, t 2 = 8. In other words, we are requiring in this example that 

£99 — £2 — t\ (52) 

With this constraint, the model can be numerically solved uniquely given the historical 
data. This assumption could be eliminated if real historical data were available for t 2 , and 
we examine the effect of varying this constraint in the sensitivity analysis. The results 
given the observed data and the constraint (52) are k a = 0.158 ± 0.066 day" 1 , k d = 1.151 ± 
0.377 day" 1 , k' = 0.122 ± 0.014 day" 1 ; the uncertainties are one standard deviation. The 
model was fit using the Frontline Systems (Incline Village, NV, USA) Solver add-in for 
Microsoft Excel 2008 for Macintosh. To obtain estimates of parameter uncertainties, we 
assumed unit variance for the input data t lt t 2) and t 99 . We then fit the sums of squared 
errors as polynomial functions of the parameters k a , kj, and k\ obtained their derivatives, 
and approximated the variances of the parameters as twice the inverse of the second 
derivative of the error with respect to each (neglecting covariances), as in [38]. We can 
now use these results as our baseline and proceed to add a hypothetical pediatric trauma 
center to this example as part of our sensitivity analysis. 

II. Sensitivity analysis 
A Approach 

In general, the output of a mathematical model depends upon the model methodology 
and the input parameters. Accordingly, the sensitivity of the output to the uncertainties 
in the parameters fit to experimental data can be assessed in a formal sensitivity analysis. 
For kinetic models of this type, much work has been published in the physics and chem- 
istry literature on methods to perform this analysis [39-42], but in this section we follow 
Atherton et al.'s approach [41]. In this section of the paper, we apply this methodology 
to fits obtained with the maximum capacity model in the previous section. In addition, 
though literature values are not available for some of the parameters in the more com- 
plicated maximum capacity with PTC model, we shall also make predictions about the 
effects of the availability of a pediatric trauma center on triage and discharge times if 
some reasonable assumptions are made about these parameters. Lastly we shall address 
mortality of the surge population using a modification of the model that includes explicit 
death rates of each population and is fully derived in Appendix B. 
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In the maximum capacity model, there are three parameters, k a , k d , and k\ and three 
outputs, ti, t 2 and t 99 . In general, if covariances are neglected, the variance or squared 
standard deviation of the /th output X b in terms of the parameters ppi a model, is: 

2 

ft (53) 



In our case, the sum index ; runs from one to three for the three parameters for 
each output variable. Therefore, there are nine elements of the relevant sensitivity 
matrix S, with the matrix elements given by 

„ _dXt 
* = 9ft 



(54) 



We can then define an output variance matrix V with the matrix elements 



(55) 



Again following Atherton et al, the effects of parameter uncertainties on the ith out- 
put variable are then ranked in order of their magnitude. 

To accomplish this, we require the nine partial derivatives implied by Equation 54, 
which are shown below: 



1 



1 
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9 £99 dt 2 



(64) 



with the first partial derivative in the numerator of the bracketed expression in Equa- 
tion 62 given by 



dk a 



3fe fl V dk a dk a 



(65) 



where the first term in Equation 65 is 



jNd ^ = ( e -k d h _ e -hh\ ( 
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(66) 
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The other partial derivative in the bracketed expression in Equation 62 is 

dN a , t2 



dk a 



fed 



(67) 



(fed - fea) 2 

Similarly, the required expressions to evaluate Equation 63 are 
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(70) 



(kd ~ kaY 



Substituting the appropriate values for the rate constants and outputs from the maxi- 
mum capacity model gives 



(-4.36 -1.14 0 
-12.23 -0.06 -49.37 
-7.87 -2.40 -49.37, 



(71) 



where the first row gives the derivatives for t v the second for t 2 , and the third for t 99 , 
and the columns correspond to differentiation with respect to k ai /c^and k\ respectively. 
After squaring each element and multiplying each column by the variance of the 
appropriate parameter, we finally obtain 



'0.082 0.184 0 
0.649 0.0005 0.500 
0.268 0.815 0.500, 



(72) 
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B. Effect of varying the constraint t 99 - t 2 

As noted in our introduction of Equation 52, in order to obtain a unique solution for 
the fit of the maximum capacity model to the data available from Reference [1], we 
had to impose a constraint on the difference between the time at which 99 percent of 
the patients had been discharged and the time at which the patient surge was 
exhausted. We arbitrarily assumed that this difference would be equal to the time 
needed to evolve from time zero to steady state, t\ in the maximum capacity model 
To determine the effect of relaxing this constraint, we varied this difference by ± 50%, 
i.e., we defined 



x = £99 — £2 



(73) 



we varied x from 1 to 3 days and re-fit the data, bracketing our initial constraint of 2 
days. The net effect of this approach is to vary t 2 , because t 99 is fixed at 10 days by the 
historical data. The results of this calculation are shown in Figure 4, which depicts the 
three rate constants k a , kj, and K as a function of x. Since t x is also constant and fixed 
by the ratio of k d to k a , as /^decreases with increasing x, /c^must increase accordingly. 
Because t 99 is fixed, and by Equation 20 the behavior of the maximum capacity model 
in region III is governed by /c^, a smaller /^results in a larger x and a shorter time 
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Figure 4 Effect of changing the constraint on t = t2-t1 on the fitted values of the rates for the 
maximum capacity model. Error bars represent one standard deviation; see text for details. 
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spent in the steady state regime of region II. Therefore, the steady state discharge rate 

K must therefore increase with increasing x, which is indeed the case. 

C. Availability of a Pediatric Trauma Center speeds admission of the pediatric cohort 

We are now in a position to include the hypothetical effect predicted by the maximum 
capacity with PTC model of the availability of a pediatric trauma center upon the flow 
of pediatric patients in this historical example. The approach we take is to vary the 
three parameters k papi k pdp , and k p from much less than the corresponding adult cen- 
ter parameters k paal k pda , an d k a f smoothly up to the latter values fitted from our his- 
torical example. We then determine the effect on observable quantities t 2 and t 99 from 
the model, the times needed to completely admit and discharge the surge population, 
respectively. For this paper, we did not independently vary the three parameters from 
zero to the fitted adult values. Rather, we first chose to look at a subset of the para- 
meter space, that in which the pediatric parameters are uniformly scaled by a single 
factor, ranging from much less than one up to nearly one, multiplied by the corre- 
sponding adult parameters. Our rationale in this approach was that without historical 
data for the ratios of the pediatric admission and discharge rates to one another, it was 
reasonable to fix them to the proportions between those of the adult center, for which, 
in contrast, we were able to fit available data. At this point, we also recall that in the 
derivation of the model, we assumed that the steady state discharge rates k p and kj 
were less than their corresponding discharge rates prior to achieving maximum capa- 
city, kpdp^id k P dai which restricts the parameter space available to explore, though this 
had no effect on the analysis that follows. 

Figure 5 shows the effect on t 2 and t 99 of varying the pediatric parameters from a 
factor of 10" 3 times the fitted adult parameters up to a factor of 0.999, and some clear 
behavior emerges. It can be seen that despite the monotonic decrease in t 2 as the 
pediatric centers effect is scaled up from near zero to approaching that of the adult 
center (Figure 5A), there is an initial increase in t 99 that peaks at a scale factor of 
approximately 0.04, and this only falls below the baseline value of 10 days when the 
pediatric parameters are scaled by 0.4 or greater (Figure 5C). We hypothesized that 
this effect arose largely from trapping of patients in the pediatric center when it was 
unable to discharge them at a sufficient rate. To test this, we investigated a second 
case where the pediatric discharge rate was fixed at the adult rate for all values of 
k paa 2Lnd k a ' t and the latter two were scaled as in the first case. As shown in Figure 5B 
and 5D, if k pdp is set equal to k pda the prolongation of t 99 is eliminated and both t 2 and 
t 99 decrease as /c^and k a ' are scaled from near zero to the adult values. 

The initial increase of t 99 for small uniform scale factors can be explained in greater 
detail by examining the behavior of the population of discharged patients in the maxi- 
mum capacity with PTC model at long times when this factor is small. In this case, we 
can write: 



where s < < 1. The population of discharged patients, the final expression in Equa- 
tion 50, can then be approximated at long times by 



kpdp - skpdt 



la 



(74) 



p d ^ i _ jj e -kpdp(t-t 2 ) 



(75) 
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Scale factor 

Figure 5 Time needed to admit (t2, panels A and B) and definitively treat to discharge (t99, panels 
C and D) the pediatric surge population for two subsets of the input parameter space. When the 
pediatric center's admission and discharge rate parameters are uniformly scaled up from zero to the values 
for the adult center, t2 decreases monotonically but t99 is initially prolonged (A, C). If the PTC discharge 
rate is set equal to that of the adult center for the pediatric surge patients while the admission and 
steady-state discharge rates are scaled up from zero to the values for the adult center, t2 (panel B) behaves 
similarly as in A, but t99 decreases uniformly (D). See text for details. 



since the C term decays much faster than the D term. It is then easy to show that 
1 , /0.01\ 

^99 = t 2 - — — In [ ) (76) 
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The derivative of t 99 with respect to s is then 
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(77) 



We must show that the right hand side of Equation 77 is positive for small but finite 
positive scale factor s. Although D is a nonlinear function of k p d p (cL Equation 43) and 
therefore of s in this approximation, its behavior is constrained by physical considera- 
tions that allow for a simple justification of this hypothesis. First, since D is the propor- 
tion of inpatients admitted to the pediatric trauma center after steady state has been 
achieved in Region III, it can never be negative, and it must necessarily be identically 
zero if the rate of admission to the PTC is also zero, or equivalently, if s vanishes. Sec- 
ondly, for very small but finite positive s < < 0.01, calculations reveal that D is positive 
but also much less than 0.01. These conditions guarantee that the second term in Equa- 
tion 77 is positive for very small s. In turn, because D increases from zero for any finite 
s, its logarithm must also increase, and the third term is also therefore positive for small 
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values of the scale factor. We note that since t 2 decreases monotonically with s (cf. 
Figure 5A), the first term in Equation 77 is negative. Despite this, numerical computa- 
tion of the values of these three terms reveals that the latter two positive terms are larger 

in magnitude than the first for small s, and dominate the behavior of — — such that t 99 

ds 

initially increases, as shown in Figure 4C. 

D. Systematic numerical sensitivity analysis of maximum capacity with PTC model 

Because we did not have historical data with which to fit the maximum capacity with 
PTC model, we chose to perform the formal sensitivity analysis assuming that the 
pediatric rates were equal to those obtained from our fit for the adult center. We set 
the variance of each pediatric parameter to 35 percent of its value, and the adult para- 
meter variances were set to the previously fitted values. We then performed the sensi- 
tivity analysis for the four outputs £ 1>0 , £ 1>/? , t 2 and t 99 as a function of the six 
parameters k paal k pap , k pda , k pdp , k a ' and k p using the same procedure as described 
above. However, for the matrix S, all partial derivatives were evaluated numerically by 
incrementing each parameter by ± 0.001, and the average value for positive and nega- 
tive increments was used for each matrix element S^. The resulting variance matrix for 
the maximum capacity with PTC model is then 

/ 0.016 0.200 0.095 0 0 0 \ 
0.016 0.200 0.095 0.217 0 0 

0.078 0.996 0.003 0.002 0.022 4.645 ( ' 

\0.027 0.345 0.888 0.907 0.023 4.690 / 

where the first row gives the magnitudes of the effects upon t lfU of changing k paa , 
kpapf k P dai k p d pi ka and k p \ the second the same values for t 1>p , the third for t 2 , and the 
fourth for t 99 . 

E Pediatric disaster-related deaths are reduced by the availability of a PTC 

A severe limitation of both the maximum capacity and maximum capacity with PTC 
models is a lack of accounting for mortality. As a final modification to the maximum 
capacity with PTC model, we included explicit death rates for each of the populations: 
the surge, pediatric patients admitted to the adult or pediatric trauma centers, and 
patients after discharge. This model is fully developed in Appendix B, and its qualitative 
behavior is demonstrated in Figure 6. We chose to use as an outcome measure the pro- 
portion of patients deceased at t = 10 days, the time at which the field hospital in Refer- 
ence 1 ceased operations. For this calculation, we began by assuming that after 
treatment and discharge, the death rate would equal the background age-adjusted death 
rate of the United States, which was approximately 8 per thousand per year in 2005 [43], 
or 2 x 10" 5 day" 1 . Although we could not find mortality data for admitted patients in the 
IDF field hospital described in Reference 1, we chose to use the figure of 8.6% mortality 
of admitted patients over 15 days, or 5.7 x 10" 3 day" 1 , from the Japanese experience after 
the 1995 Hanshin-Awaji earthquake [44]. Lastly, we based our estimate for the surge 
death rate, prior to admission and treatment, on data from the Chi-Chi earthquake in 
Taiwan in 1999, where it was reported that of all fatalities, 7% died while hospitalized 
[45] . We can therefore approximate the surge death rate by scaling our in-hospital rate 
from Reference 47 by 0.93/0.07, yielding a surge death rate of 0.076 day" 1 . We shall also 
assume that the death rate for patients admitted to the adult center is equal to that of 
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those in the PTC. With these estimates, we find that with no PTC available, the propor- 
tion of the initial surge dead at t = 10 days is 24.0 percent, but with a PTC operating 
with the same admission and discharge rates as the adult center, this is decreased to 
15.2 percent. This amounts to a reduction of the absolute mortality risk by 8.8 percent, 
and a relative mortality risk reduction of 37 percent, when a pediatric trauma center is 
available to admit and discharge patients at the same rates as those of the adult center. 

Discussion 

Summary of main results 

A deterministic first-order population kinetics model has been presented to quantita- 
tively describe the effect of the availability of a pediatric trauma center upon the time 
required to completely triage and definitively treat the pediatric cohort of a disaster 
surge. We first derived a simpler model to determine starting parameters from an histor- 
ical example. We then proceeded to examine the effect of adding in the availability of a 
pediatric trauma center over a range of values for its efficiency as described by admission 
and discharge rates relative to the baseline values obtained for the adult center. While 
the time needed to triage or admit the entire pediatric surge cohort decreased with the 
availability of a PTC regardless of its efficiency, the time to discharge of the surge had a 
more complicated behavior: if /c^is varied proportionally to the other parameters, the 
total discharge time t 99 actually increases when the PTC is slow (with rates less than 
approximately 0.04 times those of the adult center), and only begins to fall below the 
baseline value of 10 days obtained from the historical example when the pediatric rate 
constants approach 0.4 times those of the adult center. If k p d p is set equal to k p d a , how- 
ever, the times needed for admission and discharge of the entire pediatric surge cohort 
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are decreased from baseline regardless of how slow the admission or steady-state dis- 
charge rate from the PTC. Overall, if the PTC is able to admit and discharge patients at 
nearly the same rates as the adult center, the time needed to admit all pediatric patients 
is nearly halved (from 8 days to just over 4 days), and the time to complete discharge of 
the population is reduced by more than a third (from 10 days to a little more than 6.5 
days). We note that in the setting of a disaster of sufficient scale to displace a significant 
enough proportion of the population, families may not be able to receive pediatric 
patients discharged to home, so the PTC discharge rate could be retarded by this effect, 
diminishing the predicted effect on the time to discharge of the entire surge. Despite 
this, the total admission time would always be decreased with PTC availability. Lastly, 
when death rates from previous disasters reported in the literature are incorporated into 
the maximum capacity with PTC model (cf. Appendix B), we find that the overall death 
rate would be decreased from 24.0% of the initial pediatric surge population to 15.2% 
when a PTC is available to admit and discharge pediatric patients at the same rates as 
the adult center, a relative mortality risk reduction of 37%. 

The finding that t 99 initially increases when the PTC rates are uniformly scaled can be 
described as a trapping effect. In other words, when the PTC becomes available to triage 
and admit pediatric disaster surge patients, if it cannot treat and discharge them fast 
enough, then the time needed for definitive disposition of the pediatric cohort is actually 
prolonged. This occurs because overall, when the PTC is much slower than the adult cen- 
ter, the population cohort admitted to the PTC stays there much longer on average than 
those patients admitted to the faster adult center. In the context of a real disaster, this 
would result in a prolonged use of specialized pediatric hospital resources, likely increased 
costs, and a decrease in the ability of the PTC to provide routine care to the non-surge 
pediatric population. We note, however, that regardless of how slowly the pediatric surge 
cohort can be discharged, the time needed to triage and admit the surge is always 
decreased in the setting of the availability of the PTC. We therefore speculate that the 
clinical result on the surge population would be minimal, but the impairment in ability of 
the PTC to provide routine care to the background population during this period of time 
would have to be considered in disaster and contingency planning. 

Of equal interest to disaster planners are the results of the sensitivity analysis. We 
found that in the maximum capacity model (no PTC available), the discharge rate 
/c^had the greatest influence on both the time to maximum load t x and time to dis- 
charge of 99 percent of the surge population t 99 (variance matrix elements, 0.184 and 
0.815, respectively, Equation 72). In the setting where a PTC is available, the behavior 
of the total treatment time described in Figures 4B and 4D is consistent with this 
result, since the marked peaking of t 99 shown in Figure 4B is completely abolished in 
Figure 4D when the pediatric centers pre-maximum load discharge rate k p d p is set 
equal to the fitted value of /c^from the maximum capacity model in the sensitivity ana- 
lysis. On the other hand, the effect of both the pre-maximum load admission rate /c^as 
well as the steady-state discharge rate K were found to contribute about equally to the 
variance of the time needed to admit the entire cohort t 2 . These results suggest that to 
maximize the efficiency of a given center to definitively treat a given surge cohort, the 
most important factor is rapid discharge of inpatients before the maximum surge capa- 
city is reached. This observation is consistent with an analysis conducted in a large ter- 
tiary center undergoing relocation to a new facility, which found that expedited 
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discharge of inpatients was an effective means of increasing hospital capacity over the 
short term [46]. On the other hand, if the most critical goal to planners is simply to 
triage and admit the surge, with less importance placed upon definitive treatment and 
discharge, the pre-maximum load admission rate and the steady-state discharge rate 
should be optimized. 

For the maximum capacity with PTC model, the interpretation of the numerical sensi- 
tivity analysis is somewhat more complicated. For t 1}CLi the time to achieving maximum 
capacity of the faster adult center, Equation 78 suggests that the pediatric admission rate 
k pap makes the most important contribution (matrix element 0.200). This is because in 
the numerical sensitivity analysis, pediatric parameters were all varied by 35%, while the 
error in the fitted value of the adult rate k paa wdiS set equal to the much smaller error in 
kjicom the maximum capacity model. For the time to reach the maximum capacity of 
the slower pediatric trauma center, £ 1>/? , the pediatric discharge rate /c^^dominates 
(matrix element 0.217), but the pediatric admission rate /^contributes almost as much 
(matrix element 0.200). This result is not unreasonable given the explicit and implicit 
dependence of t lp upon both these rate constants (cf. Equation 43). 

In contrast, t 2 is strongly affected by the steady-state PTC discharge rate k p (matrix 
element 4.645). The dependence of t 2 on k p can be explained by two factors: first, the 
fact that the relative error in this rate assumed for our sensitivity analysis (35%) is lar- 
ger than that obtained for k a ' in the fit, and second, due to the functional form of t 2 . 
Equation 48 shows that t 2 is a linear function of t lp and B, and therefore depends 
implicitly on A^and k p a p . It is inversely proportional to the sum of the adult and 
pediatric steady-state discharge rates, k = k p + k a \ We have observed that invariably, 
whenever t Xp increases, regardless of whether k pc i a or k pc i p is changed, B decreases. 
Therefore, the effect of changing either k pda ov /coupon t 2 is limited because of this 
antagonistic effect. In contrast, changing k by varying k p or k a ' does not produce a 
compensatory change in either t Xp or B, so the effect of k p dominates. 

Lastly, t 99 depends most strongly on k p , with the next strongest dependence on A^and 
/^(matrix elements 4.690, 0.907, 0.888 respectively). Though we cannot write down an 
analytic expression for t 99 in the maximum capacity with PTC model, we can make quali- 
tative arguments based on the behavior of this parameter in the simpler maximum capa- 
city model. Equation 51 reveals that in the simpler model, t 99 depends explicitly on t 2 and 
the discharge rate k d , with implicit dependence upon both /<^and /^within the argument of 
the logarithm. It is reasonable to conclude that in the more complicated maximum capa- 
city with PTC model, the dependence would be similar on t 2 and the two discharge rate 
constants /c^and k p d a > Since we have seen that for the maximum capacity with PTC 
model, t 2 is most sensitive to changes in k p \ it follows by this reasoning that k p will also 
have a large effect on t 99 . Moreover, we have already seen that in the simpler model, /^ac- 
tually has the greatest effect on t 99 , so taken together, this combined with the qualitative 
argument discussed here provide a reasonable explanation for the sensitivity of t 99 to 
kpdpand k pda in the maximum capacity with PTC model. 

Limitations of the model 

The potential methodological weaknesses of the model must also be considered. First, 
as noted above, no distinction is made in the discharged populations of either the max- 
imum capacity model, or the maximum capacity with PTC model, between patients 
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discharged home or otherwise dispositioned, including discharge to nursing care cen- 
ters, rehabilitation hospitals, and even death. Similarly, the death of patients in the 
surge population prior to triage and admission is not accounted for. This concern is 
addressed in full detail in Appendix B, where a more complicated version of the model 
including death rates is derived, and is implemented in additional files 1 and 2. We 
expect that a deterministic population kinetics approach will describe the behavior of 
the populations of interest only when they are sufficiently large. However, for very 
small populations the continuous mathematics used to derive our model would be 
expected to break down, and a discrete stochastic approach [47] might be more 
appropriate. 

Tradeoffs 

An important consideration for disaster planners is the potential cost of various 
approaches to preparedness. Though our model provides a mathematical justification 
for the inclusion or use of a pediatric trauma center in the response to a disaster, it does 
not consider the monetary cost of establishing one, or the resources required to keep it 
in operation. The average cost of building a new hospital has been reported in the Uni- 
ted States to be approximately 285 dollars per square foot as of 2003 [48], or 342 2011 
dollars per square foot [49]. At our own facility, a new 460,000 square foot (42,700 m 2 ) 
specialized children's hospital with 317 beds, a level I pediatric trauma center and sup- 
porting facilities cost 636 million dollars in 2011, a cost of nearly 1400 dollars per square 
foot [50] . Therefore, prior to committing to building such a facility, a careful accounting 
of the likelihood of various types of disaster occurring in the proposed construction area 
as well as the availability of rapid transportation to and capacities of already existing 
nearby centers would have to be performed. Alternatively, a different approach would be 
for planners at an established center to prepare mobile dedicated pediatric trauma cen- 
ter facilities similar to the mobile field hospital described in reference 1, available to be 
transported to the site of a disaster as needed. However, we speculate that this method, 
though much less expensive than building a new PTC, could possibly have detrimental 
effects on treatment of affected adult patients. For example, after prolonged operation, 
such facility would require resupply, and if a medical resupply shipment had to be par- 
celled out to the PTC in addition to competing adult centers in the affected area, the 
resulting relative shortage of resources in the adult centers might result in decreased 
rates of admission and discharge, and increased death rates, of adults. Such considera- 
tions, though beyond the scope of our model directly, would also have to be examined 
to allow for its use in disaster planning. 

Conclusions 

The model presented here provides an analytical, closed-form description of the popu- 
lation dynamics of a disaster surge population treated either in the presence or the 
absence of a pediatric trauma center, is mathematically elementary and is simple to 
implement. Given that the proportion of children in the population is roughly twenty- 
five percent, 35 the potential influence of the availability of a specialized trauma center 
whose resources are devoted to the pediatric surge cohort must be taken into consid- 
eration by public health agencies. We have demonstrated how the model can be 
applied to an historical example to obtain starting parameters, and the hypothetical 
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contribution of an available PTC can then be assessed as a function of how it com- 
pares in efficiency to the historical example. If detailed quantitative historical data that 
explicitly included a PTC as part of the response to a disaster became available, the 
model could be fit to these data and estimates of the model parameters could be 
obtained. While the costs of building and maintaining a PTC and the effects of its 
resource consumption on other hospitals must be taken into account, this determinis- 
tic kinetic model provides a new weapon in the armamentarium of disaster planners. 
Our approach can be used to provide a hypothetical estimate of how the response to 
an historical event could have been improved, as well as to extrapolate and predict 
potential responses to future events. 

Appendix A 

General case of Eqs. 1-4 with surge delayed from inciting event 

For all MCEs, there is a delay between the inciting event or exposure and the develop- 
ment of the associated patient surge. The approximation made in the treatment in this 
paper is that the delay is much smaller than any of the other timescales in the model (i. 
e., admission or discharge). This is an excellent approximation for sudden MCEs such as 
bombings, earthquakes or airplane crashes. However, for some classes of MCE, such as 
disease pandemics, radiation exposure events, floods, hurricanes, as well as the aftermath 
of more sudden types of insults considered above, the delay time between event and 
surge is of the same order of magnitude as these other timescales, and must be treated 
explicitly in the model. 

We now present the general case only for Region I of the simpler model entailed by Eqs. 
1-4, because the inclusion of the pediatric trauma center makes the equations significantly 
more complicated, with the introduction of an additional parameter (the delay time) and 
differential equations, with a limited contribution to any further physical or planning 
insight. The procedure for obtaining the solution in the maximum capacity and zero 
queue-length regimes (Regions II and III of the maximum capacity model), as well as the 
inclusion of a PTC, would be the same as that in the main text. In the general case of the 
model, in the absence of the pediatric trauma center, we would have four populations 
rather than the three in Eq. 1: 



Here /c 5 is the exposure or delay rate; the remaining rate constants are identical to 
those of Eq. 1. Also, instead of N 0 instantaneous surge patients, we now have N 0 
exposed patients at time zero. The governing differential equations are then: 
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The boundary conditions are: 
N e (t = 0) = N 0 

N s (t = 0) = N a {t = 0) = N d (t = 0) = 0 
N 5 (t -> oo) = N fl (t -> oo) = N d (t -> oo) = 0 
The solution to (A2-A5) is then: 



N a (t) = N 0 
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Similarly to the case of Eqs. 31 and 32, we can compute exactly the time of maxi- 
mum expected surge by taking the derivative of A10 and setting it equal to zero, 
which gives: 

1 



K - h 



In 



(A13) 



We note that in the context of disaster planning, A13 can either be used to predict 
the time of maximum surge, if estimates for /c^and /c 5 are known, or to constrain and 
relate k a to k s if the maximum surge time is known from historical or data or other pre- 
dictive methods. 

Appendix B 

Maximum capacity with PTC model and explicit death rates 

In this section we shall derive a version of the maximum capacity with PTC model 
where a background death rate of each population is included. For times at which 
neither the adult nor the pediatric center has reached maximum capacity (analogous to 
region I of Figure 3) the governing differential equations are: 
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where co s , co a , <d p , and co d are the death rates for the surge, pediatric patients admitted to 
the adult center, pediatric patients admitted to the PTC, and discharged patients, respec- 
tively, Po{t) is the total number of deaths that have occurred at time £, and all other para- 
meters are as defined in the main text. We introduce /, K and L as the boundary 
conditions for the deceased population P D (t) at t h a1 t h pi and t 2 , respectively. We define 

X s = k + co s 

K = kpda + 0J a (B2) 
hp = kpdp + M p 

The solution to (Bl) is: 
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The time at which the adult center's maximum capacity is reached is 
h,a = In I — 



(B4) 



For times at which the faster adult center has reached its maximum capacity, but the 
PTC has not (analogous to region II of Figure 3), the system's behavior is governed by 
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where <D a ' is the constant death rate of patients hospitalized in the adult center, and 
we have introduced 
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The solution to (B4) is 

P s (t) = a'e-^-* 1 *) - p 

Paa (f) = C 

P ap (t) = (H + <p-6)e- x p( t - t ^ 

+Q e -hl(t-h, a ) _ y 

P d (t) = r (e~ x p^~ h ^ - g-^(t-ti,a)^ 

+ A fg-M'-^a) - e -Vd{t-h,a)^ 

+A (l - e-^-' 1 ^ + Ee-^-^ 

P D (t) = ((Da + AoJ d - pCO s - (pC0p) (t t ha ) 

+ _L ((H + <p — 0) co p + Tco d ) (l - e-^-^ 

+ {a' (o s + 0co p + AcD d ) (l - e~ kn ^~ h ^ 
+ (A + A + T - E) U-^(t-h,a) 

Where 



A = e~ kstha 
k 



C 



H 



■paa /-X s t ha _ e -k a t ha \ 



Xa — X 



k 



-pap 



Xp — X. 



1 / kp da kpaa kpdp\ 



co d - X s \X a - X s X p -X s ) y 



kpdak 



■■pan 



(C0 d - X a ) (X s - X a ) 

kpdpkpap 
(co d - X p ) (X s - X p ) 



^0—^ah,a 



-<*>dh,a\ 



0 



/=i(l - A) 



kpaa kpap 
00 s + CO a — + (On 



"Xa ~ 



Xt> — Xc 



1 / kpdafepaa kpdpkpap 
+(D d — I — + 



fe 



'Xa Q^s ~ 'Xa) 



k 



hp (X s — Xp} 

+ (l - e - mh ' a ) 

kpdpkpap 



(D d — X$ \ Xa — X$ Xp — X$ J 

(co a + (D d ^-) (1-e-^) 
V (Dd- K) 

- ) (1 -e"Vl.a) 



(Dp + (Dd 



kpdp 
(Dd ~ \ 



kpda)^) 



■pan 



1 



1 



_X S — X a \co d — X s (D d — X a 
1 1 



X s — Xp \(D d — X s cod — X t 



(B6) 



(B7) 



(B8) 



Barthel et al. Theoretical Biology and Medical Modelling 201 1, 8:38 
http://www.tbiomed.eom/content/8/1/38 



Page 28 of 32 



and 
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We note that we have introduced a new boundary condition / for the fraction of 
patients who are deceased at t h a . The time at which the pediatric center reaches its 
maximum capacity is again obtained by maximizing P ap (t) in B6 and is 
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After this time, and until the surge is exhausted, both the adult and pediatric centers 
are at steady state and can only admit a patient if another is discharged or dies: 
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and (Dp' is the constant death rate of patients hospitalized in the pediatric center dur- 
ing this time period. The solution to B9 is 
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where K is the fraction of deceased patients at t h p and the constants are given by 

B = a'e~ hl ^~ h ^ - p 
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The time t 2 at which the surge is exhausted is determined by setting the first expres- 
sion in Bll equal to zero and is 

t 2 = h p + — lnf 1 +B— ) (B15) 

For times after t 2 the adult and pediatric centers resume discharges at pre-maximum 
capacity rates: 

dP s 

t>t 2 : — - = P s (t) = 0 
at 

— - —X P 

dt " a aa 



dP ~ -XpPap (B16) 



dp, 

= + Vp^ap - COdPd 

dP D 

= 0) aa Paa + MpPap + COdPd 

The solution to B13 is 

p aa (t) = ar^ (£ - £2) 

P ap (t) = De^-^ 

p d (t) = C _^_ fe-A*(t-fe) _ ^(t-fe)) 

m - K 

+D _^pdp_ < -X p (t-t 2 ) _ e -co d (t-t 2 )\ 
COd - Xp 

+C e -vd(t-t 2 ) (B17) 



Xp\ Q)d~ XpJ 

(1 _ e -K(t-t 2 )^ 



C ( CQdkpdc 
+ — COa + 



( G C kpdg _ D kpdp \ 
V CDd-ka COd-XpJ 



e -cod(t-t 2 )\ 



Barthel et al. Theoretical Biology and Medical Modelling 201 1, 8:38 
http://www.tbiomed.eom/content/8/1/38 



Page 30 of 32 



where L is the fraction of deceased patients at t 2 and the constants not yet defined 



Additional material 



Additional file 1: We shall briefly describe here the first Microsoft Excel spreadsheet that accompanies 
this paper. The spreadsheet, titled MC and MCPTC models.xls [Additional File 1], allows the reader to enter values 
for the rates of admission and discharge, both for the maximum capacity model (worksheet tab "Max capacity 
model") and for the maximum capacity with pediatric trauma center model (worksheet tab "Max capacity model 
with PTC"). The behavior of each population and resulting timescales in the two models, as well as the boundary 
conditions described in the text and illustrated in Figure 3 for the latter model, are calculated and displayed in the 
labelled boxes as described by header notes on each worksheet.Regarding the calculations themselves, the usual 
Excel worksheet functions EXP and LN are employed, but in order to properly calculate the behavior of the 
populations in each regime, conditional logic statements must be constructed. In computer programming 
languages, statements of the form iLthen; else... must be coded in such cases. In the case of the maximum 
capacity model, for example, this involves evaluating different expressions for each population for the time 
regimes r <fi, fi < f <f 2 , and f > t 2 . To accomplish this in Excel, the worksheet function \F{logical test, expression 7, 
expression 2) can be constructed to contain nested conditions [51]. In other words, since we require a statement 
of the formif t <t ± then f x (f) / else if f x <, t <t 2 then f 2 (f) ;else f 3 (t) ;end if where f^t), etc. are the 
functions that describe the desired behavior in each region, the corresponding Excel statement islF("f <t 1 ", "fi(0", IF 
(AND("r > U", "t <t 2 "), "f 2 {t)", "f 3 (r)")where the quotation marks reflect the fact that the enclosed statements are just 
shorthand for demonstration purposes, and in an Excel worksheet must be properly formatted Excel statements 
and functions of worksheet cells. For a demonstration please see the additional files.We also reiterate here that the 
derivation for the maximum capacity with PTC model (section III of Methods in the manuscript) assumes that the 
pediatric trauma center is, at most, no faster than the adult center. This constrains the values that can be input by 
the user for the discharge rate of pediatric patients from the PTC to be less than or equal to that of their 
discharge from the adult center, i.e. k pdp < k pda . If the user sets k pdp >k pda , the admitted population of the pediatric 
center reaches its maximum capacity before the adult center, and the derivation (which assumed the opposite is 
true) is invalid. To explore the behavior of the model when the PTC is faster (not addressed in this paper), the user 
could set the pediatric rate constants to the desired or known baseline adult center values, and then vary the 
adult rates to be as rapid as desired. This amounts to simply relabeling all the rate constants for the PTC as those 
for the adult center and vice versa. 

Additional file 2: The second Excel spreadsheet [Additional File 2] incorporates both the maximum 
capacity with PTC model and explicit death rates for each pediatric population as derived in Appendix B. 

The user interface is similar to that of the first file, with modifiable inputs in bold and boxes color-coded for inputs 
(blue), dimensionless constants (green), and outputs and boundary conditions (orange). The lower right hand 
corner of the orange box also displays the proportions of discharged and deceased patients. Because of the 
greater complexity of this version of the model, we also included worksheet tabs for regions I, II, III, and IV 
separately for demonstration purposes, which can be found to the right of the main model tabs labelled "MCPTC 
with deaths, PTC unavail" and "MCPTC with deaths." 
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